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ABSTRACT 

We present the mass functions for different mass estimators for a range of cosmological models. We pay 
particular attention to how universal the mass function is, and how it depends on the cosmology, halo 
identification and mass estimator chosen. We investigate quantitatively how well we can relate observed 
masses to theoretical mass functions. 

Subject headings: cosmology: theory - large-scale structure of Universe 



1. INTRODUCTION 

One of the most fundamental predictions of a theory of 
structure formation is the number density of objects of a 
given mass, the mass function. Accurate mass functions 
are used in a number of areas in cosmology; in studies 
of galaxy formation, in measures of volumes (e.g. galaxy 
lensing) and in attempts to infer the normalization of the 
power spectrum, the statistics of the initial density field, 
the density parameter or the equation-of-state of the dark 
energy from the abundance of rich clusters. One of the 
most intriguing aspects of the mass function is that it ap- 
pears universal, in suitably scaled units, for a wide range 
of theories. A complete understanding of this phenomenon 
currently eludes us. 

If we are to attempt to use the mass function to infer 
cosmological parameters from the abundance of objects 
of some given property, then we need to understand how 
accurate our theory for the mass function is. This in- 
volves understanding how to define the mass of an object 
in cosmology, a task which is non-trivial as there is no 
clear boundary between a halo and the surrounding large- 
scale structure in theories of hierarchical structure forma- 
tion. The purpose of this paper is to calculate the mass 
functions for different mass estimators for a range of cos- 
mological models, to see how well these statistics can be 
computed from semi-analytic theories and to investigate 
quantitatively how to relate observed masses to theoreti- 
cal mass functions. 

We find that the mass function is only approximately 
'universal', and only for a few mass estimators. We discuss 
how accurately one can convert between different mass es- 
timates, so as to relate what can easily be measured to 
what can easily be predicted. We also discuss the limita- 
tions which non-universality of the mass function would 
place on cosmological parameter estimation were it not to 
be corrected for. 

The outline is as follows: after a review of some back- 
ground (§^j) we present mass functions, derived from N- 
body simulations (§|J), for a variety of different mass def- 
initions (§|). We compare these mass definitions, based 
on mean density contrast, with the concept of a virializcd 
halo (§||). We investigate how universal the mass function 
is (§^) for each different estimator and present fitting func- 
tions (§Q) to the mass functions for 3 different cosmological 
models. We finish by considering the effect of clustering 
on the mass function (qq) and summarize our main results 
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2. PRESS-SCHECHTER THEORY 



We begin by reviewing the basic theory underlying the 
expectation that the number density of halos, per unit 
comoving volume, should take a universal form. Thi s ex- 



pectation was first elucidated by Press & Schechter f974 
who combined the statistics of the initial density field with 
a model for the evolution of perturbations based on spheri- 



cal collapse of a top-hat overdensi ty (see e.g. Peacock |I999 
for a textbook treatment; Bower |199f|; Peacock & Heav- 
ens 



1990 



Bond et al. 1991 Lacey & Cole 1993 for more 
details). Specifically these authors advanced the ansatz 
that the fraction of mass in halos more massive than M is 
related to the fraction of the volume in which the smoothed 
initial density field is above some threshold 5 C . A variety of 
smoothing windows and thresholds have been advocated, 
but the most common is a top-hat window in real space 
and 5 C ~ 1.69. 

The P-S mass function agrees relatively well with the 
results of numerical simulations both for critical density 
models with power-law spectra and, more surprisingly, for 
model s with out this self-similar evolut ion (e.g. Efstathiou 
et al. [1988|; E fstathiou & Re es |1988l; White, Efstathiou 
& Frenk |1993|; Lacey fc Co le |1994| ; Celb & Bertschinger 



1994| ; BoncTS Myers [1996D - tneP-S mass function and 



numerical results are known to deviate in detail at both 
the high and low mass ends. Refinements to this theory 
have been advanced, all of which relate the abundance of 
collapsed objects to peaks in the initial density field in 
a 'universal' manner. In the latest incarnation the mass 
function has been motivated by or fit to l arge cosmological 
N-body simulations (Sheth & Tormen 1999; Jenkins et 
al. |200C| ; hereafter JFWCCECY). 

To fix our notation we recap briefly the ingredients in 
this model in the next two sections. 

2.1. Top-hat collapse 

The spherical top-hat ansatz describes the formation of 
a collapsed object by solving for the evolution of a sphere 
of uniform overdensity S in a smooth background of den- 
sity p. By Birkhoff 's theorem the overdense region evolves 
as a positively curved Friedman universe whose expansion 
rate is initially matched to that of the background. The 
overdensity at first expands but, because it is overdense, 
the expansion slows (relative to the background) and even- 
tually halts before the region begins to recollapse. Tech- 
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nically the collapse proceeds to a singularity but it is as- 
sumed in a "real" object virialization occurs at twicefj] the 
turn-around time, resulting in a sphere of half the turn- 
around radius. In an Einstein-de Sitter model the over- 
density (relative to the critical density) at virialization is 
A c = I871- 2 ~ 178. We shall always use A c to indicate the 
overdensity relative to critical of a virialized halo, which 
will be lower for smaller Cl m . A fitting function for A c 
for arbitr ary f2 m and £Ia can be found in Pierpaoli, Scott 



& White 2001[ Note that some authors use a different 
convention in which A c is specified relative to the back- 
ground matter density - our A c is f2 m times theirs and we 
shall come back to this point in §|4]. The linear theory ex- 
trapolation of this overdensity is normally denoted S c and 
is (3/20)(12tt) 2 / 3 ~ 1.686 in an Einstein-de Sitter model. 
This overdensity is often used as a threshold parameter in 
PS theory and its extensions and has a very weak cosmol- 
ogy dependence which is often neglected. We shall return 
to some of these considerations in ^ 

2.2. Multiplicity function 

The mass function now comes from considering the 
statistics of the initial density field and the top-hat model 
above. Under the P-S ansatz all of the cosmology de- 
pendence is contained within the rms density fluctua- 
tion, er(M), smoothed with a top- hat filter on a scaleQ 
R 3 = 3M/47rp. The multiplicity function, 



„. , , Af dn 

vf{v)dv = dM 

Jy ' 2pdM 



(1) 



is a universal function of the peak height v which is related 
to the mass of the halo through 



(2) 



a(M) 



with 5 C = 1.69. Not e that some authors, particularly 
Sheth & Tormen |1999| , define v to be {S c /a(M)) 2 rather 
than 6 c /<j(M) as we have done. If the initial fluctuations 
are Gaussian, as we shall assume throughout, then the 
multiplicity function is simply 



vf(v) oc e 



(3) 



where the normalization constant is fixed by the require- 
ment that all of the mass lie in a given halo 



J vf{v)dv = ~ 



(4) 



There is no justification for this normalization from N- 
body simulations, which cannot probe the M — > tail, 
but we shall adopt it throughout. 

M otivat ed by a model of elliptical collapse, Sheth & Tor- 
men 1999 provided a fit to large, high-resolution N-body 



simulations of the modified form 



f{v) oc (1 + (ais 2 )-?) (av 2 )-V 2 e- a » 2 / 2 



(5) 



1 There is a small correction to this in the presence of a cosmolog- 
ical constant which contributes a Ar 2 potential. 

2 In principle R could be defined with respect to p cr i t , but this is 
not the natural choice in the top-hat collapse model. 



where p = 0.3 and a — 0.707 provided the best fit to 
groups selected with a spherical overdensity algorithm. A 
slightly different fitting function was proposed by JFWC- 
CECY based on analysis of the same simulations. The 
advantage of Eq. (§) over that of JFWCCECY is that it 
is well behaved over the full range of mass, whereas the 
functional form of JFWCCECY cannot be safely extrap- 
olated outside of the range of their fit. In addition the 
elliptical collapse model can be used to discuss the clus- 
tering of halos (Sheth & Tormen 1999) using an extension 
of the p eak-backgroun d spli t formalism ( e.g. E fstathiou et 
al. |198S| ; Cole & Kaiser |l989j ; Mo & White |1996| ). While we 
shall see later that the Sheth & Tormen form overpredicts 
the number of small mass halos in some of our simulations, 
that region is not the main focus of this work. Small dif- 
ferences in the functional form of the fitting functions will 
not be important for our conclusions. 

2.3. Toward higher accuracy 

To recap the material in the last 2 sub-sections, we as- 
sume that the mass function of virialized objects depends 
only on the variance of the initial density field, smoothed 
on some scale with a specified filter. We calculate the frac- 
tion of the volume which is occupied by peaks which exceed 
a threshold value and relate this to the number density of 
halos of a specified mass. The critical density threshold is 
taken from the theory of spherical top-hat collapse and is 
treated as a constant. 

At first sight it seems fortuitous that any result of the 
complex process of halo formation within an hierarchical 
model (see e.g. Fig. |l|) could be derived from the variance 
of the initial density field, without re ferenc e to any dynam- 
ics (see discussion in Lacey & Cole 1994 ). Or that there 
should be a 'universal' mass function at all. Indeed JFWC- 
CECY found that the cosmology- independence of the mass 
function in scaled units depends upon the mass estima- 
tor chosen. This is a non-trivial problem because objects 
in hierarchical models do not have a well defined outer 
boundary, making both their identification and the defi- 
nition of their total mass convention dependent. JFWC- 
CECY found best results when using as a mass estimator 
the sum of the particle masses in their N-body g roups us- 
ing a particular group finder (FoF; Davis et al. 1985, see 
§0). This differs from the widely followed practice of us- 
ing the mass within a spherical region w hos e radius i s de- 
rived from the top-hat collapse model (§2.1). White 2001 



showed that there is considerable scatter between the two 
types of mass estimators which makes it difficult to com- 
bine results which don't use a consistent set. 

Is there a middle ground? The self-similarity of halos 
observed in simulations has long been taken to imply that 
masses arc best defined within radii enclosing a fixed den- 
sity contrast. While the density contrast A c is the con- 
ventional choice, it is not the only possibility. Since FoF 
links particles which are approximately above some den- 
sity threshold with respect to the background, the result 
of JFWCCECY suggests that we should define our masses 
within fixed density contrasts with respect to the mean 
density, not the critical density. JFWCCECY in fact give 
such a mass function in their Appendix B. It uses a mass 
within a radius rrsob interior to which the mean density is 
180 times the background density. While this is close to 
the 'top-hat' result for an O m = 1 cosmology, it extends to 
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extremely large radius compared to the observable region 
of clusters. For example r 180 b ~ 2 — 3/i _1 Mpc for a rich 
cluster (M - 10 15 h^ 1 M Q ). 

A number of cosmological tests rely on the existence of a 
mass function which is both universal and easy to interpret 
observationally. None of the results described above obvi- 
ously fulfill these two requirements. We shall try to make 
steps towards this goal in the rest of this paper. Our final 
solution will be a hybrid which uses the mass estimator 
Misob suggested by JFWCCECY along with a conversion 
factor between 'observed' and 'theoretical' mass. 

3. SIMULATIONS 

Numerical simulations give qualitative support to the 
predictions of the Prcss-Schcchter theory, with small mod- 
ifications noted at both the high and low mass ends. The 
current state of the art in numerical simulations aimed at 
elucidating these departures is the work of JFWCCECY. 
Independent confirmations of the J FWC CECY result have 
been published recen tly by White [200]] , Zheng et al. |2002 
and Hu & Kravtsov 2002, In this section we discuss the 
numerical simulations we have done to investigate the de- 
pendence of the result on the mass definition chosen. The 
reader not interested in the numerical details is urged to 
skip to §f| 

3.1. N-body runs 

We have performed a suite of N-body simulations in 
order to better constrain the mass/multiplicity function 
(see Appendix for details of the code). The first set we 
used to tune our fitting function, the rest were used as 
independent checks. Throughout we have tried to focus 
primarily on the high-mass end of the mass function, which 
is of the most use for studies of clusters of galaxies. 

Since the primary consideration is one of volume, we 
have run numerous small simulations rather than one very 
large one. The small simulations were chosen to have suf- 
ficient dynamic range and mass resolution to well resolve 
a low mass cluster halo. Specifically we ran a number of 
150 3 particle simulations of three different CDM models, 
each in a 200/i _1 Mpc box (see Table pi for more details). 
Each simulation represents a reasonable cosmological vol- 
ume, so as not to bias the high mass end of the mass func- 
tionrL while maintaining enough mass and force resolution 
to identify the relevant halos. 

Because it provides a reasonable fit to a wide range of 
observations, we first simulated a 'concordance' ACDM 
model which has O m = 0.3, fi A = 0.7, H = 100 /ikms _1 Mpc~ 
with h = 0.67, n B = 0.04, n = 1 and cr 8 = 0.9 (corre- 
sponding to Sh = 5.02 x 10 -5 ). We call this model 1. 
We then changed the model in each of two 'orthogonal di- 
rections'. First we changed the mapping between length 
scale and mass, by changing il m from 0.3 to 1, while hold- 
ing the present day power spectrum fixed (model 2). Then 
we changed the normalization of the power spectrum, erg , 
while holding the cosmology and shape of the power spec- 
trum fixed (model 3). Finally we ran a model with a dif- 
ferent spectral shape and normalization as a cross check. 
Model 4 had n m = 0.35 = 1 - fi A , V B h 2 = 0.02 and 
h = 0.65 with eg = 0.8. Of the models the first has the 

3 The contribution from modes with wavelength longer than the 
box to cr(R) in the relevant range is very small. 



fewest clusters per unit volume, so we ran more realiza- 
tions of this model than the other three. 

We also used an 'independent' simulation to check our 
fitting function. The cosmology in this case was slightly 
different than above, having f2 m = 0.3, h — 0.7 and as = 1. 
The simulation employed 512 3 particles in a 300ft- _1 Mpc 
with a smoothing length of 20/i _1 kpc. We can use this 
simulation to see whether the mass function extrapolates 
correctly to lower masses (where it becomes a power-law) , 
see §|. 

3.2. The group catalogs 

From the z — output of each simulation we produce 
a halo catalogue by ru nning a "friends-of-friends" group 
finder (e.g. Davis et al. 1985| ) with a linking length of ei- 
ther b = 0.2 (in units of the mean interparticle spacing) 
or b = 0.1. We can use these two different group cata- 
logs to test the sensitivity of our results to the selection 
of halos. The FoF algorithm partitions the particles into 
equivalence classes, by linking together all particle pairs 
separated by less than a distance b. We keep all groups 
above 32 particles, which imposes a minimum halo mass 
of order lO 13 /i -1 Af . [FoF groups with more than 32 par- 
ticles are known to be robust.] The FoF algorithm as we 
have defined it cannot be used to address sub-structure in 
the halos that we find, but for our purposes this will not be 
a serious limitation as the P-S formalism also completely 
neglects sub-structure. 

Several other group finders exist which we could have 
used in addition to the FoF algorithm. Some of these be- 
gin with groups defined in a FoF manner while some find 
a partition of the particles in a completely different way. 
Luckily, an exploration of all of these different group find- 
ers will turn out to be unneccesary. For most of the mass 
estimates defined in §^ the precise group finding algorithm 
is unimportant. We shall show later that the mass func- 
tions obtained with two different FoF groups, with radi- 
cally different partitionings of the particle distribution, are 
almost identical. This may be telling us that the physical 
properties of the group which we are calculating, the 'to- 
tal mass', are independent of the details of how the group 
is originally found as an overdensity in three dimensional 
space. 



Model n r 



o~8 N hc 



frpart 



0.30 
1.00 
0.30 
0.35 



0.9 
0.9 
1.0 
0.8 



15 
10 
10 
10 



120 
80 
80 
80 



1.97 
6.58 
1.97 
2.30 



40 
30 
40 
40 



Table 1 

The parameters of the simulations run. In each case 
nbox different realizations of the gaussian initial 
conditions were run. each run used a periodic box 
of side 200fo _1 mpc and 150 3 particles. the force 
softening was of a spline form with a "plummer 
equivalent" smoothing length of 50/l _1 kpc - easily 

SMALL ENOUGH TO RESOLVE THE HALOS OF INTEREST. In 
THE ABOVE VOLUMES ARE QUOTED IN UNITS OF 
(100 /l -1 Mpc) 3 AND PARTICLES MASSES IN UNITS OF 

lO u ft _1 M . 
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Fig. 2. — The projected density in a cube 10 h~ x Mpc on a side centered on the second most massive halo in the 512 3 particle simulation. 
The 3 panels are projections down the x, y and z axes of the box. The grey scale is logarithmic, running from 10 2 to 10 s times the mean 
density. The solid circles show r200c — 1.74/i -1 Mpc (inner) and rigot = »*54c — 3.04h -1 Mpc (outer). Within rigot the material exhibits a 
wide range of density contrasts. Note that the halo is neither isolated nor spherical, and has quite a bit of substructure. 
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In order to define the mass it will be very useful to have 
a halo center. We define the center of a halo as the position 
of the potential minimum, calculating the potential using 
only the particles in the FoF group. This proved to be 
more robust than using the center of mass, as the poten- 
tial minimum coincided closely with the density maximum 
for all but the most disturbed clusters. Additionally, the 
position of the center was very insensitive to the presence 
or absence of the particles near the outskirts of the halo, 
and thus to the precise linking length used. 

4. THE MASS OF A HALO 

As remarked earlier, since the objects formed in a hier- 
archical model have no clear boundary, all mass definitions 
are a matter of convention. For each halo in each catalog 
we computed 8 different definitions of mass. As we shall 
see later, it is only those definitions which encompassed 
the majority of the virialized material ([Q) which will turn 
out to be useful for our purposes, and the best shall be 

^1806- 

First we used the 'FoF mass', simply the number of par- 
ticles in the group times the particle mass. With b = 0.2, 
this definition is the one used by JFWCCECY, for which 
they found a universal mass function. By considering the 
mean number of particles in a sphere of radius b one can 
argue that (if all particles have the same mass) FoF groups 
are bounded by a surface of density 3/(27rfe 3 ) ~ 60 times 
the background. If all groups were spherical and singular 
isothermal spheres^], this would imply a mean density in- 
side the FoF group of roughly 180 times the background 
density or 180f£ m times the critical density. In practice 
there is a very large scatter about this value. We also use 
the sum of the particles in the b = 0.1 groups. 

Motivated by the self-similarity exhibited by halos in 
simulations we also define the mass from a spherically av- 
eraged profile about the cluster center. Specifically we de- 
fine Ma as the mass contained within a radius r& inside 
of which the mean interior density is A times the critical 
density 

J r 2 dr p(r) = — p C rit?A . (6) 

The 'virial mass' from the spherical top-hat collapse model 
would then be simply Ma c . We shall refer to this mass as 
Mth-vir rather than the somewhat vague term virial mass. 
Since A c ~ 200 in a critical matter density cosmology, 
many authors mean by 'virial mass', M^oo- This is in fact 
the most common definition. We shall write this M200C to 
make explicit the fact that it is with respect to the critical 
density. We shall also consider M5000 which has the ad- 
vantage that it probes material at sufficiently small radii 
that it is often directly accessible to X-ray observations. 

In the above we have followed common usage and mea- 
sured the mean interior density contrast to the critical den- 
sity. This has historically been motivated by (a) consid- 
erations based on the virial theorem, which provides esti- 
mates of the halo velocity dispersion and 'temperature', in 
which the critical density provides a natural scale and (b) 
because it requires no assumptions about the cosmological 
parameters (beyond h) in its definition. More recent work, 
specifically JFWCCECY, has suggested that the halo mass 

4 If p(r) oc r~ 2 the mean density interior to a radius where the 
density is p is just 3p. 



function may be universal if masses are measured within 
a fixed density contrast measured not with respect to the 
critical density, but with respect to the background den- 
sity. We shall follow their lead and also calculate rigob, 
where the mean interior density is 180 times the back- 
ground density. Note that while the mass function may 
be more universal with this definition, it comes with an 
associated price from an observational point of view: it 
introduces a dependence on an assumed f2 m in the defini- 
tion. 

Unfortunately both the FoF mass with b = 0.2 and risob 
encompass a very wide range of material (see Fig ||), far 
beyond what can usually be observed and into the region 
where the profiles start to show significant scatter from 
halo to halo. For a cluster mass halo, r^ob can be 75% 
larger than the most widely used r2ooc- For this reason we 
shall also consider M§oob and Miooofc, which require less of 
an extrapolation. 

In cases where we use a small linking length and a large 
ta it is possible that two halos overlap. Since our defini- 
tion of mass for each halo includes all of the mass within 
ta, not just that associated with the FoF halo itself, this 
can result in us double-counting the mass in the overlap 
region. This is an unfortunate side-effect of a mass esti- 
mator based on spherical averages for objects which are 
neither spherical nor isolated. To avoid this we cull from 
each run the smaller of two halos whose centers are closer 
than the sum of the 'virial' radii. This procedure is then 
relatively insensitive to the linking length used to define 
the original halos. If we had used a larger linking length 
and 'merged' the extra halo with the larger one (removing 
it from out list), the mass assigned to the combined halo 
would still be that within ta of the potential minimum of 
the combined group - presumably the potential minimum 
of the larger mass group - and thus be the mass originally 
assigned to the larger of the two groups. If we do not 
perform this culling step then we find a significant excess 
of low mass objects compared to the analytic predictions. 
With the culling the mass function becomes almost totally 
insensitive to the original group finder parameters. 

The need for this step is more than just a technical is- 
sue. It stems from the fact that in hierarchical models 
halos are rarely isolated, but are often found in various 
stages of merging or accretion. Observers often remove 
from their samples systems which they deem to be inter- 
acting too strongly or too recently. This can introduce a 
'bias' in the mass function. Rather than attempt to quan- 
tify how 'disturbed' various halos are or whether there is 
observational evidence for interaction which would cause 
them to be removed from any particular sample, we have 
chosen to apply a criterion that can be equally well ap- 
plied in simulations and observations: we omit the smaller 
of any two systems whose virial radii overlap. The frac- 
tion of the halos culled in each of the models is relatively 
small. For example in Model 1 for the b = 0.2 halos the 
culled fraction is less than 1% for all of the mass defini- 
tions. For b = 0.1 it is 15% for Mi 80 &, 8% for M th - V ir 
and around 1% or less for the other mass definitions. How 
closely our treatment mimics the selection of individual 
objects in current samples is unclear. As we demand in- 
creasing accuracy in our comparison between theory and 
observation this issue will need to be revisited. 
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5. THE VIRIAL REGION 

Our definition of a halo is primarily one of density con- 
trast. An alternative definition is that a halo contains ma- 
terial which has broken away from the universal expansion 
and is (at least approximately) in virial equilibrium. As 
we shall sec below, it is only for density thresholds where 
these two definitions roughly coincide that one obtains a 
nearly universal mass function. Since the virial region en- 
compasses such a large volume of space, this will require 
us to investigate extrapolations from the inner, easier to 
measure, regions. 

As a first step it is interesting to ask how well the spheri- 
cal top-hat collapse model approximates the messy forma- 
tion process of a cluster in a hierarchical model (see Fig. Q). 
To investigate this we have looked at the five most mas- 
sive clusters in the 512 3 particle simulation described ear- 
lier. For each cluster we calculate the mean radial velocity, 
v r (r), and the velocity dispersion, <r(r), in bins containing 
5000 particles from the center out to 3 times the virial ra- 
dius predicted by the top-hat model (rioi for Q m = 0.3). 
The results are shown in Fig. ||[ 

Fig. H shows that the clusters are roughly isothermal, 
with a velocity dispersion profile that peaks near the break 
radius (where the density profile has a slope of —2). In- 
side of the virial radius the mean velocity is close to zero, 
in units of the 3D velocity dispersion er vir . Just outside 
the virial radius the transition from 'virialized' material 
to in-flowing material (v r < 0) is clearly seen in all 5 clus- 
ters. Thus it seems that the 'virial radius' is predicted 
within a factor of 2 by the top-hat model, though the 
clusters exhibit some scatter. The radius ri 80 b is only 30% 
larger than r t h-vir for a cluster mass halo in this cosmol- 
ogy, so within the scatter we could take the 'virial' radius 
to be risoh also. Using a higher density contrast than A c , 
e.g. M200C or M5000 would clearly underestimate the tran- 
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Fig. 3. — The radial velocity (dotted) and velocity dispersion 
(dashed) profiles of the 5 most massive clusters in the 512 3 simu- 
lation described in the text. Each cluster has 0(1O 5 ) particles and 
-M20OC > 10 15 h~ 1 Mq. Points are plotted every 5000 particles in 
radius out to 3x the top-hat virial radius. The velocities are all 
normalized to the 3D velocity dispersion of the dark matter within 

r th-vir- 



sition radius for all of the clusters. 

To get a feel for the translation between density con- 
trast and 'size' for a rich cluster, we can make use of the 



universal density profile of Navarro, Frenk & White 1996 
These authors defined the virial radius as r2ooe arid den- 
sity contrasts with respect to critical. The universal form 
of the density profile is 



f 



p{r) 

oc — — 

Pent X(l+Xy 



(7) 



where x = r/r s is a scaled radius and r s describes the tran- 



sition from 



to r 



in the profile. We show the radius 



within which the mean density is A times the critical den- 
sity in Fig. [| for a rich cluster with M 2 oo c = 10 15 h~ x M Q 
and a concentration c = r200c/r s = 5. Note that for 
Sl m = 0.3, risob = ^54 C ~ 2.8 h~ 1 Mpc. For lower fl m it 
would exceed 3/i -1 Mpc. 

6. A UNIVERSAL MASS FUNCTION? 

Given a (possibly culled) set of halos each with a 
known mass, we wish to find a fitting function to the 
mass/multiplicity function. As a first step we check 
whether the mass functions do indeed form a 'universal' 
multiplicity function by plotting 



N(> v) 



M dn 
p dv 



dv 



(8) 



vs. the peak height v 2 (Fig. ||). As we can see, for the 3 
models showrj^J, it is a good approximation to assume that 
all of these mass functions come from a universal multi- 
plicity function for the FoF halos with b — 0.2. The uni- 
versal form is quite well fit by the Sheth & Tormen form 
of Eq. (§). This confirms the earlier work of JFWCCECY. 

The mass function is also close to universal if the top-hat 
virial mass is used, though the agreement is not as good 



5 We omit here the results of Model 4 for visual clarity. 
Model will be reinstated in some later plots. 
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1000 



Fig. 4. — The radius r& within which the mean density is A times 
the critical density for an NFW halo with M20OC = 10 15 h~ 1 MQ and 



'T>4, 



5. For f2 m = 0.3 the 
2.8/i~ 1 Mpc. 



universal' density contrast of r 1 g 0b is at 



7 




Fig. 5. — The multiplicity function vs. the peak height v 2 for our 3 models and for 6 of our 8 mass definitions. Open circles are Model 1, 
open squares Model 2 and triangles Model 3. The solid line is the S-T fitting function. The panels are labeled with the halo mass definition 
used. This indicates to what extent the mass functions arc indeed universal. 



8 



as in the b = 0.2 case and would need to be checked for a 
wider range of cosmologies. As noted by JFWCCECY the 
mass -M18O6 also gives a close to universal mass function. 
Each of these mass estimators includes the majority of the 
mass within the virialized region (see Fig. |^). 

As expected, the mass function within r2ooc shows a 
systematic difference between the J7 m = 1 model and the 
other two. This is because the mass is defined interior to a 
density contrast which doesn't scale with the background 
density. 

Of particular interest is the case of Miooofc- Here the 
mass is defined in terms of a density contrast with respect 
to the background density (like the 'universal' Migofc), 
but at a higher density. This region of the halo is more 
amenable to observation. However we note that as we in- 
crease the density threshold, focusing on the inner regions 
of the halos, the universality of the mass function degrades. 
[The case Msoob is intermediate, and is omitted from the 
figure.] 

We can understand this result by considering the dif- 
ferent formation histories of the halos. Fig. ^ shows the 
average mass profiles of the most massive halos from the 3 
models in scaled units. The halos in Model 2, with fi m = 1, 
form later and thus are less concentrated than the halos 
in Models 1 and 3. This increases the ratio Miooofc/-^i80fc 
(see Fig. ffj) or Mi oob/M 2 oo c - On an object by object ba- 
sis this scatter is fairly large. This is unfortunate since 
it is precisely the inner regions which are most amenable 
to observation! Thus the mass estimators for which the 
mass function is close to universal are those which re- 
quire an extrapolation beyond the observed region (out 
to 2 — 3/i _1 Mpc), and this extrapolation is quite sensi- 
tive to both the cosmology and the particular formation 
history of the halo under consideration. 

We show in Fig. || the mass functions on a linear scale. 
We report the results as a ratio of the N-body results to 




Fig. 6. — The (interior) mass profiles in scaled units for the halos 
above Migot = 5 X 10 14 /i -1 Mq in the 3 models. As before, open 
circles are Model 1, squares Model 2 and triangles Model 3. The 
mean and standard deviation of the profiles is shown. Note that by 
definition M(< riggb) = 180M(< r^gob) so there is no scatter in 
this point. 
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1 I000b/^180b 



Fig. 7. — The ratio of different mass estimators for clusters with 
M 18 ob > 5 X f0 14 h^ 1 Mq in the 3 models. Solid line is Model 1, 
dotted Model 2 and dashed Model 3. 
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Fig. 8. — The multiplicity function vs. the peak height u 2 7 divided 
by the fitting function of Eq. (a) , for Mxgot. ■ Open circles and stars 
are Model 1 with b = 0.2 ando = 0.1 respectively. Open squares 
and 3-pointed crosses Model 2, triangles Model 3 and plusses and 
diamonds Model 4 (with worse statistics at the high mass end). 
The 4-pointed crosses are the results from the 512 3 run of a slightly 
different model (see text). At the high mass end (v — > 00) our 
statistics become very poor. 
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the integral of Eq. (H) for the case of Migob- Each of 
the four models is represented by two sets of points, one 
where the halos are initially found with a linking length of 
b = 0.2, and the other with b = 0.1. As we can see there is 
almost no dependence on the initial group finder used in 
this case. As a cross check we also show on this linear scale 
the mass function from the 512 3 particle simulation in the 
300/i _1 Mpc box. This run has much higher mass and 
force resolution than the majority of the runs used here, 
and a slightly different cosmology i.e. spectral shape and 
normalization on the scales of interest. The total volume 
is however smaller, so the statistics at the v — + oo end 
are much poorer. However it makes a good 'independent' 
check of the universality of the mass function. The scatter 
between the models is at the ±20% level. 

We also show, on the same linear scale, the results 
of converting between different sp herica lly averaged mass 
profiles in Fig. ^. Following White 2001 we show in partic- 
ular what happens if one constructs the mass function by 
measuring M 500c , converting this to Mi 80 b using an NFW 
profile with c = 5. We make no correction for the large 
scatter we saw in the detailed comparisons above, we sim- 
ply apply a numerical rescaling. For Q ro = 0.3 the conver- 
sion is -M18O6 — 2.OA/500C while for Q m = 1 the conversion 
is Migob — lAMsooc. As we can see, while the conver- 
sion shows a lot of scatter it introduces no major bias and 
the mass functions so constructed are close to universal. 
Such a pro cedure was follo wed, in reverse, by Pierpaoli 
et al. 2001 and Seljak 2002 for example. A more compli- 
cated conversion a long the same lines has been provided 
by Hu & Kravtsov ^002| . This is very encouraging because 
it means that, while Misofc cannot reliably be measured 
on an object-by-object basis, a noisy estimator of it can 
be easily constructed which turns out to be good enough 
to construct the 'universal' mass function. The outlier 
is Model 4 which was already known to be discrepant in 
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Fig. 9. — The multiplicity function vs. the peak height u 2 , divided 
by the fitting function of Eq. (ph as above. Here we have converted 
from Msooc to MigQi, assuming the halos are all of the NFW form 
with c = 5. Open circles and stars are Model 1 with b = 0.2 and 
b = 0.1 respectively. Open squares and crosses Model 2, triangles 
Model 3 and plusses and diamonds Model 4. 



Fig. |[ The rescaling has made it more discrepant from 
the mean, indicating that this procedure is not without its 
flaws, but even so the mass function is predicted to 30% 
over much of the range. 

7. FITTING THE MASS FUNCTION 

For completeness we would like to find a fit to the sim- 
ulated mass functions (Figs. [l0|, [ll]). Since to a very good 
approximation the mass function is independent of the 
clustering of the halos, we can do this using the Poisson 
model. First we bin the halos in mass using a large enough 
number of bins that no bin contains more than 1 halo. We 
use bins equally spaced in logM. Then we maximize the 
(log) likelihood 



log£ = lo S^i 



iefuii 



fij + const 

j Gall 



(9) 



where the sum on i is over bins containing 1 particle, the 
sum on j is over all the bins and 



dn 



dlogM 



A log M < 1 



(10) 



is the mean number of halos per bin assuming a mass 
function dn/ dlogM. This method has the advantage of 
being independent of the chosen bin width and correctly 
taking into account the Poisson statistics of the rare halos 
at the high mass end. 

We have chosen to use the modification to the Press- 
Schcchter formula, Eq. (||), proposed by Sheth & Tor- 
but to allow a and p to be free parameters. 



1999 



We maximize the likelihood for a and p fitting over the 
range M = 10 14 h- l M Q to M x = 3 x lO^h^Mg. The 
results of this procedure are shown in Table |[ We found 
that the fit was somewhat sensitive to the particular value 
of Mq chosen, indicating that the numerical mass func- 
tion wasn't perfectly fit by the form of Eq. (||). Using 
the higher resolution simulation we found that the fitting 
function tended to overpredict the abundance of halos with 
v 2 < 1/2, by a factor of almost 2 for the lowest mass halos 
we could probe. Since we concentrate here on the more 
massive end of the mass function we did not attempt to 
correct for this. 

There is a fair degree of variation in the best fit pa- 
rameter a, while p is close to constant. This is because 
p essentially controls the slope of the low-mass (y ~ 1) 
end of the mass function which is very nearly the same for 
all the estimator s. We obtain reasonable agreement with 
Sheth & Tormcn 199£ for the conventional mass estimator 



M200C for the critical density cosmology for example, but a 
is significantly smaller using the b = 0.2 FoF mass and sig- 
nificantly larger using the b = 0.1 FoF mass. Generally a 
increases as the mass estimator probes material primarily 
at a higher density contrast. 

For the 'universal' contrast of Migob calculated directly 
from the halo mass distribution we found better agree- 
ment with the numerical results in all cases if we lowered 
a slightly below the 0.7 found by Sheth & Tormen for this 
mass estimator. This is the opposite of the claim of Hu 



& Kravtsov 2002 that a should be increased to 0.75 when 
using Mi806- On the other hand for the one example we 
studied in detail where we converted from a mass mea- 
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sured within r$ooc to one within rigoh using an NFW pro- 
file with c = 5, the best fitting a was slightly higher than 
0.7. The difference in the mass functions so produced, 
as a is changed from 0.65 to 0.75, is at approximately 
the same level of the scatter from model to model shown 
in Fig. ||. Thus these fluctuations could be reflecting an 
intrinsic limit to how well we can determine our fitting 
function. 

8. CLUSTERING AND THE MASS FUNCTION 

Throughout we have assumed that the number of clus- 
ters in a given volume is simply Poisson distributed about 
a mean value given by the mass function. In principle 
however the positions of clusters are correlated an d this 



can induce non-Poisson fluctuations (Evrard e t al. 2002 



for a simple analytic model see Hu & Kravtsov 2002 ) . We 
expect this effect to be small when the objects are rare and 
the sample region is large compared to the cluster corre- 
lation length r ~ O(10 /i^Mpc) (see e.g. Peebles |1980| ) 
but we can quantif y its e ffect using numerical simulations. 

Hu & Kravtsov 2002 have investigated the additional 
scatter in the normalization that arises from including 
clustering using a simple analytic model where clusters 
are biased tracers of the linear density field. We shall in- 
vestigate this effect using the z = group catalog from 
the Hubble Volume Simulation of a ACDM model run by 
the Virgo Supercomputing Consortium (JFWCCECY). To 
quantify the scatter in the mass function, we shall com- 
pute the best fitting power spectrum normalization, erg, 
to 1000 random sub-volumes of the simulation, using the 
maximum likelihood method described above. We hold 
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Table 2 

Mass function parameters for the different 
cosmologies and different mass definitions (see 
text). 



the other parameters fixed at their fiducial values for sim- 
plicity. Each sub-volume is centered on a random point 
in the simulation, which can be chosen as the center us- 
ing the periodicity of the box. Then all clusters are kept 
which have M above a threshold mass, are closer than R 
to the center of the box and would have |6| > 30° if the 
box x — y plane was oriented parallel to the galactic plane. 
This roughly mimics a volume limited X-ray survey out to 
depth R sensitive to the most massive, and therefore most 
clustered but rarest galaxy clusters. We choose two mass 
thresholds, 10 14 h^ 1 Mq and 3 x W 14 h' 1 M G . The scat- 
ter expected from a Poisson distribution can be estimated 
using the same procedure, except that we first randomize 
the positions of the halos. 

Fig. |lj shows the standard deviation in erg, in units of 
the mean value, as a function of sampled volume for the 
'clustered' and 'Poisson' cases. We checked that estimat- 
ing the variance directly or from the difference between 
the 16 th and 84 th percentiles of the distribution gave sim- 
ilar results. For all volumes studied the variance in as is 
increased by cluste ring over the simple Poisson expecta- 
tion (Evrard et al. 2002), however for volumes of interest 
(-Rmax > 300 /i^MpcJboth the Poisson variance and the 
increase due to clustering are almost negligible co mpare d 
to the other errors (see e.g. Table 6 of Pierpaoli et al. 2001). 
We checked explicitly that there was no bias in the mean 
introduced by the neglect of clustering in the analysis. 

9. CONCLUSIONS 

The multiplicity function, a measure of the number of 
halos per comoving volume element per unit mass, is one 
of the central predictions of a model of structure forma- 
tion. Dark matter dominated models in which structure 
evolves hierarchically from gaussian initial conditions pre- 
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Fig. 12. — The standard deviation of erg, in units of the mean, as 
a function of survey radius. We plot a quantity, half the difference 
between the 16 th and 84 th percentiles of the distribution, which is 
slightly less sensitive to outliers than the variance; but the results 
would be almost identical if we had plotted the variance. Solid lines 
indicate the width of the distribution including the clustering of 
clusters, dashed lines are for the randomized sample. Lines joining 
open squares are for M > 3x 10 14 Mq and joining solid triangles 
for M > lO^h- 1 M P) . 
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Fig. 10. — Mass functions for the 3 models, 2 linking lengths and 4 of the mass definitions. Left panels are for b = 0.2, right for b = 0.1. 
Top to bottom are Models 1, 2 and 3. The open squares are FoF mass, solid triangles M th _ vil , open circles M20OC and crosses Msooc- In the 
first panel the solid line is the Press-Schechter prediction, and the dashed line the fit to the Virgo simulations. The solid line is suppressed in 
the other panels. 
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Fig. 11. — Mass functions for the 3 models, 2 linking lengths and 4 of the mass definitions. Left panels are for b = 0.2, right for b = 0.1. 
Top to bottom are Models 1, 2 and 3. The open squares are FoF mass, solid triangles M lgoi) , open circles M 5 Q Qb and crosses M 1000b . In the 
first panel the solid line is the Press-Schechter prediction, and the dashed line the fit to the Virgo simulations. The solid line is suppressed in 
the other panels. 
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diet a mass function which is nearly universal if expressed 
in the right units. This is only true for a narrow class 
of mass estimators, and is specifically not true for the es- 
timators which have been most commonly used up until 
now. 

Given the complicated process by which halos form in 
hierarchical models, the role of mergers and prevalence of 
sub-structure, it is highly convenient that the mass func- 
tion (in scaled units) is so close to universal. We have pro- 
vided fitting functions to the mass function from N-body 
simulations for 8 different mass estimators, and shown how 
one can convert between them. We have found that mea- 
suring the mass of a halo using one definition and using a 
simple average spherical profile, such as the NFW profile, 
to convert to the 'universal' Migob provides a remarkably 
good method of estimating the mass function, even though 
individual halos show a large scatter among different mass 
estimates. 

Finally let us remark upon the small non-universality in 
the multiplicity function, which can lead to misestimates 
of the true mass function if one uses a fitting function 
like Eq. (||). Neglecting the factor (ilog cr/cilog M in the 
mass function, making an error of Sn/n in the number 
density per logM at mass M translates into an error of 
(y 2 — l)~ 1 8n/n in the normalization da/a. So for a typical 
cluster, with v ~ 2 — 3, the scatter in the mass function 
doesn't limit our knowledge of og until the other uncertain- 
ties are pushed below 0(5%). Thus the non-universality of 
the mass function is not currently a limitation to using the 
abundance of rich clusters to determine the normalization 
of the matter power spectrum. The uncertainties become 
increasingly important when it comes to using the evolu- 
tion of the mass function as a probe of f2 m or the equation 
of state, w, of the dark energy. For the latter, errors on 
ag approaching the percent level are required. For these 
ambitious measurements it may not be sufficient to use 
a simple parameterized form for the multiplicity function. 
One could either resort to full blown numerical simulations 
for a grid of models 'near' the parameter region of inter- 
est or attempt to find a 'second variable' which correlates 
well with the scatter between the simulation results and 
the P-S predictions. As we approach the level of preci- 
sion where these effects matter a variety of other effects 
also become important, i nclud ing the effects of clust ering 
(see e.g. §|, Evrard et al. ^002| , Hu & Kravtsov |2002j ) and 
how to treat merging systems. It will be a challenge for 
theorists to keep the 'theory uncertainty' below the 'exper- 
imental uncertainty' with the increasingly rapid advances 
in observations. 
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THE TREEPM-SPH CODE 

The simulations in this paper were done using the 
TreePM-SPH code (White et al. |2002| ) running in fully 
collisionless mode. We present a brief discussion of the 
features of the code here for completeness. 

The TreePM code was specifically designed to run on 
distributed memory computers or clusters of networked 
workstations and evolves dark matter (and gas) in a pe- 
riodic simulation volum e. The code uses the TreePM 
method of Bagla |1999 for the gravitational force and 
smoothed particle hyd rodyn amics (SPH; Lucy 1977 and 
Gingold & Monaghan 1977 ) in its 'entropy' formulation 
(Springel & Hernquist [2002 ) to compute the hydrody- 
namic forces. The (collisionless) dark matter component 
and (collisional) gas are assumed to interact only through 
gravity. The code is written in standard C and uses the 
Message Passing Interface (MPI) communications pack- 
age, making it easily portable to a variety of parallel com- 
puting platforms. It performs dynamical load balancing 
and scales efficiently with increasing numbers of proces- 
sors. It will run on an arbitrary number of processors, 
though it is slightly more efficient if the number is a power 
of two. 

The particles are integrated using a second order leap- 
frog method, where the relevant positions, energies etc are 
predicted at a half time step and used to calculate the 
accelerations which modify the velocities. The time step 
is dynamically chosen as a small fraction (depending on 
the smoothing length) of the local free-fall time. Particles 
have individual time steps so that the code can handle a 
wide range in densities efficiently. To increase speed the 
force on any given particle is computed in two stages. The 
long-range component of the force is computed using the 
PM method, while the short range component is computed 
from a global tree. In this manner the code is similar 
in spirit to P 3 M except that the short range force, being 
computed from a tree, scales as NlogN rather than N 2 . 

Rather than a Plummer potential we use a spline soft- 
ened force. We use the same smoothing kernel for the 
gravity and SPH calculations. The long-range force is 
smoothed on 2 grid cells, and the opening criterion for the 
tree is set to achieve 1% accuracy in the short-range force. 
With these standard parameters the 90th percentile force 
error is 1.2% for lightly clustered distributions while for 
very uniform distributions the 90th percentile error rises 
slightly to 1.9%, as the total force is smaller and the short 
range force contributes less. 

We have made extensive comparisons of the code de- 
scribed here to other codes described in the literature, 
building on the many test problems that those codes have 
been shown to satisfy. In particular we have compared 
TreeP M ex tensively with Gadget (Springel, Yoshida & 
White 2001). Details of these comparisons, along with re- 
sults from self-similar evolution, hydrodynamics tests in- 
cluding the sod shock tube, th e San ta-Barabara cluster 
comparison project (Frenk et al. 1999] ) etc can be found in 
White, Springel & Hernquist |2002j. 



The Virgo data are publicly available at http://www.mpa- 
garchng.mpg.dc/NumCos. This work was supported in 
part by the Alfred P. Sloan Foundation 
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